Linear and nonlinear superparamagnetic relaxation 
at high anisotropy barriers 
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The micromagnetic Fokker-Planck equation is solved for a uniaxial particle in the low-temperature 
limit. Asymptotic series in the parameter that is the inverse barrier height-to-temperature ratio are 
derived. With the aid of these series, the expressions for the superparamagnetic relaxation time and 
the odd-order dynamic susceptibilities are presented. The obtained formulas are both quite compact 
and practically exact in the low (with respect to FMR) frequency range that is proved by comparison 
with the numerically-exact solution of the micromagnetic equation. The susceptibilities formulas 
contain angular dependencies that allow to consider textured as well as randomly oriented particle 
assemblies. Our results advance the previous two-level model for nonlinear superparamagnetic 
relaxation. 
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I. INTRODUCTION 

The problem of superparamagnetic relaxation in 
single-domain ferroparticles formulated, explained and 
basically analyzed by Neel[[l] about fifty years ago, since 
then keeps to attract attention. Nowadays this inter- 
est is mainly due to the expanding number of nanometer 
granular magnetic media used in information storage and 
related high technologies. 

When analyzing magnetic dispersions, solid or fluid, 
a promising idea is to evaluate the granulometric con- 
tent, particle material parameters and relaxation rates by 
combining the data on linear and nonlinear dynamic sus- 
ceptibilities. Since recently, this approach (it originates 
from the spin glass science) became quite feasible in ex- 
perimental realization.]^] However, to benefit from it, one 
needs an adequate model. Surprisingly, until nowadays 
the NeelQ concept of superparamagnetic behavior of fine 
magnetic particles that had been substantially advanced 
by Brown || ^| and refined by numerous researchers (see 
the review article [|| with about 400 references), lacks a 
nonlinear extension. 

In Ref. |^ we begun to fill up this gap and proposed 
a numerical procedure involving continuous fractions by 
means of which the linear and cubic susceptibilities for a 
solid system of uniaxial fine particles could be obtained. 
With allowance for the polydispersity of real samples, the 
worked out description provided a fairly good agreement 
with the dynamic magnetic measurements taken on Co- 
Cu nanocomposites.pl Recently, our approach was used 
successfully for the linear and cubic susceptibilities of 
the samples of randomly oriented 7-Fe203 nanoparticlcs. 
Hereby we carry on the build-up of the nonlinear super- 
paramagnetic relaxation theory by working out a set of 
compact and accurate analytical expressions that consid- 
erably facilitate calculations as well as experiment inter- 



pretation. 

The paper is arranged in the following way. In Sec. H 
we discuss the problem of superparamagnetic relaxation 
and show the way to obtain the asymptotic solution for 
the micromagnet ic F okker-Planck equation in the uniax- 
ial case. In Sec. Ill the perturbative expansions for the 
orientational distribution function are obtained, which 
are used in Sec. IV to construct asymptotic expressions 
for the nonlinear dynamic susceptibilities. The explicit 
forms of those expansions are given and their accuracy 
is proved by comparison with the results of numerical 
calculations. Sec. ^ contains the enveloping discussion. 



II. SUPERPARAMAGNETIC RELAXATION 
TIMES 

A. Uniaxially anisotropic particle 

The cornerstone of the superparamagnetic relaxation 
theory is the Arrhenius-like law for the relaxation rate 
of a magnetic moment of a single domain particle pre- 
dicted by Neel in 1949. The framework of this classical 
problem is as follows. Consider an immobile (e.g. fixed 
inside a solid matrix) single-domain grain of a volume 
v. This particle possesses a uniaxial volume magnetic 
anisotropy, K being its energy density and n its easy 
axis direction. Since the temperature T is assumed to be 
much lower than the Curie point, the particle magneti- 
zation /, as a specific parameter, is practically constant 
and the magnitude of the particle magnetic moment may 
be written as fi — Iv. Denoting its direction by a unit 
vector e, one concludes that the magnetic state of such 
a particle is exhaustively characterized by a pair of vec- 
tors: n = Ive and n. Thence, the orientation-dependent 
part of the particle energy (in the absence of external 
magnetic fields) is 



U = -Kv(e-ri) 2 , 



(1) 
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where K is assumed to be positive. As Eq. (|]J) shows, 
this energy has two equal minima. They are separated 



2 



by the potential barrier of the height Kv and correspond 
to e || ±n because for the magnetic moment e the direc- 
tions n and — n are equivalent. At zero temperature, the 
magnetic moment e, once located in a particular poten- 
tial well, is confined there forever. At finite temperature, 
the probability of an overbarrier (interwell) transition be- 
comes non-zero. If the ratio a = Kv/kT is high enough, 
the transition rate is exponential thus yielding the Neel 
law r tx exp(cr) for the reference time r of the particle 
re- magnetization. 

BrownQ shaped up those semi-qualitative considera- 
tions into a rigorous Sturm-Liouvillc eigenvalue problem 
by deriving the micromagnetic kinetic equation 



2r D dW/dt = JWJ(U/kT + lnW), 



(2) 



where W(e,t) is the orientational distribution function 
of the magnetic moment, J = (ex d/de) is the infinites- 
imal rotation operator with respect to e, and the time 
td is introduced below by formula (^). Generally speak- 
ing, equation (|^) is incomplete since a gyromagnetic term 
is absent there. This means that the consideration is 
limited by the frequency range lutq <C 1, where To is 
the relaxation time of the Larmor precession of the par- 
ticle magnetic moment in the internal anisotropy field 
H a ~= 2K/I, where K includes the possible shape con- 
tribution. Comparing this condition with the other one, 
^lTq ^$ 1, that evidences a low-to- moderate quality fac- 
tor of the Larmor precession for real nanodisperse fer- 
rites, one estimates the allowed frequency as lo <C lol 
that means, in fact, a fairly wide range. ]25| 

In the statistical description delivered by Eq. (g), the 
observed (macroscopic) magnetic moment per particle is 
given by the average 



m(t) = n(e) = J eW(e,t)de 7 



Note that with allowance for Eq. ([y) the function W has 
a parametric dependency on the vector n so that, in fact, 
the angular argument of W is (e • n). 

The magnetodynamic equation underlying the Brown 
kinetic equation ^ can be either that by Landau & Lif- 
shitz or that by Gilbert. To be specific, we adopt the for- 
mer one. Thence, the reference relaxation time in Eq. (|^) 
writes as 



td = Iv/2ajkT, 



(4) 



(3) 



where 7 is the gyromagnetic ratio for electrons and a 
is the precession damping (spin-lattice relaxation) phe- 
nomenological parameter. 

Assuming uniaxial symmetry of the time-dependent 
solution and separating the variables in Eq. (||) in the 
form 

W(e,t) = ± J2Zo ^ Me • n) exp(-A^/2r rJ ), (5) 

where the amplitudes Ag depend on the initial perturba- 
tion, one arrives at the spectral problem 

Ltf^i = Ae4>t, L = J 2cr(e ■ n)(e x n) — J . (6) 
where the non-negativity of the decrements can be 
proven easily. Expanding the cigenmodes tpi m the Leg- 
endre polynomial series 



1> t = \ £(2fc + l)&f P fc (cos0), k= 1, 3, 5, ... , (7) 
fc=i 



where 9 is the angle between e and n, one arrives at the 
homogeneous tridiagonal recurrence relation 



A/ 



k(k + l) 



2a 



k-1 



(2k- l)(2fc + l) 



'k-1 



1 



(2k- l)(2fc + 3) 
I 



(2/s + l)(2fc + 3) 



J k+2 



(8) 



Note that Eqs. (||)-(||) describe only the longitudinal 
(with respect to the easy axis) relaxation of the magnetic 
moment. We remark that under condition lu <C oj-l, i.e., 
far from the ferromagnetic resonance range, the transver- 
sal components of m = /i( e ) are of minor importance. 



B. Interwell mode 

Spectral equation (^|) describes the temperature- 
induced (fluctuation) motions of the vector e in the ori- 
entational potential with a symmetrical profile ([!]). With 
respect to the time dependence, the set of possible eigen- 



modes splits into two categories: interwell (overbarrier) 
transitions and intrawell wanderings. In the spectral 
problem (^) the interwell transitions of the magnetic mo- 
ment are associated with the single eigenvalue Ai. As the 
rigorous analysis shows, Jj| it drastically differs from the 
others: whereas for £ > 1 all the \g gradually grow with 
a, the decrement Ai exponentially falls down proportion- 
ally to exp(— a). 

In the opposite limit a — >■ 0, all the decrements, in- 
cluding Ai, tend to the sequence \p, — i(i + 1) and thus 
become of the same order of magnitude. This regime 
corresponds to a vanishing anisotropy so that the differ- 
ence between the inter- and intrawell motions disappear, 
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and the magnetic moment diffuses almost freely over all 
the 47r radians with the reference time tjj introduced by 

From Eqs. (Bh and (H) one finds that the longitudinal 
component of the magnetic moment evolves according to 



i(t) = pb Aj e 



-\ e t/2T D 



1=1 



xipi dx, 



(9) 



where a; = cos (9 = (e ■ n). For a symmetrical potential 
like the equilibrium value mo of the particle magnetic 
moment is zero. 

With the above-mentioned structure of the eigenvalue 
spectrum, the term with i = 1 in Eq. ([}]), being propor- 
tional to exp(— e~ a t/T£,), at a > 1 is far more long-living 
than any other one. The dominating role of the decre- 
ment Ai had been proven by Brown, and for it he had 
derived 0] the asymptotic expression 



A B = (4/Vtt) o 



3/2 



(a » 1). 



(10) 



In a short time after, using a continued fraction method, 
Aharoni constructed || for Ai a fairly long power series 
in a and also showed numerically that Brown's expres- 
sion ([n]) resembles the exact one with the accuracy of 
several percent for a > 3. In the 90 's the eigenvalue 
Ai became a subject of extensive studies. Efficient nu- 
merical procedures were developed and a number of 
extrapolation formulas with a good overall accuracy were 
proposed || (i| 0. 



C. Asymptotic solution of the Brown equation 

The study that we describe below was inspired by 
our work on fitting the dynamic susceptibilities mea- 
surements for real assemblies of fine particles. Those 
data typically describe polydisperse systems in the low- 
frequency bandwidth oj/2it — l-10 3 Hz. As tq ~ 10~ 9 s 
or smaller, then, using formula ( |l0| ) for estimations, one 
concludes that the mentioned frequency interval becomes 
a dispersion range for the interwell (superparamagnetic) 
mode at 



LUTqC 



> 1 that is a> 10. 



For temperatures up to 300 K this condition holds for 
quite a number of nanomagnetic systems. 

Application of the best fit procedure to a set of ex- 
perimental data implies numerous re-calculations of the 
linear and non-linear susceptibility curves x^ of the as ~ 
scmbly. Any such curve, due to a considerable polydis- 
persity of the particles, is a superposition of a great num- 
ber of partial curves X (°0 spread over a wide size (or, 
in the dimcnsionlcss form, a) range. For successful pro- 
cessing, one needs a fast and very accurate algorithm to 
evaluate x^ (<?") everywhere including the domain a 3> 1 . 
The existing extrapolation formulas are no good for that 
purpose due to their ill-controllable error accumulation. 
A plausible way out is an asymptotic in ct -1 solution of 
Eq. (JgJ) . In the course of the fitting procedure, this ap- 
proximation can be easily matched in the intermediate 
a range with the well-known expansions for the small a 
end. 

It is noteworthy that some 20 years ago Brown him- 
self resumed[[l5| |l6| studies on Ai and modified the pre- 
exponential factor in Eq. ( |l0| ) transforming it into an 
asymptotic series in a^ 1 . On the base of Eq. @ he had 
constructed an integral recurrence procedure, and eval- 
uated Ai down to terms cx 1/er 10 . What we do below, 
is, in fact, carrying on this line of analysis that had not 
been touched since then. Our method advances Brown's 
results in two aspects. First, for Ai it is more simple. 
Second, it provides not only the eigenvalue but the eigen- 
function as well. Only having the latter in possession, one 
is able to obtain theoretical expressions for the directly 
measurable quantities that is the susceptibilities x^ ■ 

Taking Eq. (^) as the starting point, we remark its 
equilibrium solution 



= Zq 1 exp(ra 2 ), Z = 2R(a), (11) 

-^( cr ) = / exp(cra; 2 ) dx, 
Jo 



that corresponds to I = and Ao = and note the 
asymptotic expansion for the partition integral R(cr) 
found in Ref. fi~7l 



R(*) = e°G/2*, G(.).l + ^ + A + ^ 

I 



(2n - 1)!! 

2 n a n 



(12) 



The operator L in Eq. (||) is not self-conjugated and 
thus produces two sets of eigenfunctions, which obey the 



respective equations 



Ltpk = AfeVfc 



(13) 
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here + denotes Hermitian conjugation. The eigcnfunc- 
tions of these two families are orthonormalized and re- 
lated to each other in a simple way: 



dxifjipk = 5jk 



(14) 



Qualitatively, from Eq. (|l4|) one may say that (p). are the 
same eigenfunctions but "stripped" of the exponential 
equilibrium solution ip . Substituting Eq. (|l4|) in Eq. (|^), 
one gets two useful relationships 



ipo(J<Pj)(J<fk) dx = XkSjk , 
(15) 

where the second one follows from the first after multipli- 
cation by ipi and integration by parts. Note that in the 
second formula action of each operator reaches no farther 
than the nearest closing parenthesis. 

On rewriting Eq. (|l5|.l) in terms of a single orienta- 
tional variable x = (e • n), the spectral problem takes the 
form 



_d_ 

dx 



V>o(l-* ) 



2\ dtpk 



dx 



-Xkipatpk- 



In the equilibrium state Eq. ( |ilf ) reduces to 

= 0, 



d 

dx 



2 ^ dipo 



(16) 



(17) 



whose normalized solution is c/'o = 1- This solution, be- 
ing a true equilibrium one, turns the inner part of the 
brackets, i.e., the probability flux in the kinetic equation 
(0), into identical zero. 



As remarked in Sec. II B, at a 3> 1 the most long-living 
non-stationary solution of Eq. ( |l6|) is the eigenfunction 
with I = 1, whose ei gen value is exponentially small, see 
Brown's estimation (|l0|). We use this circumstance for 
approximate evaluation of ipi in the (7>1 limit by ne- 
glecting the right-hand side of Eq. ( jlrj ) for 1 = 1. On 
doing that, the equation obtained for the function ipi 
formally coincides with equation ( p~T| ) for (fQ. However, 
the essential difference is that now the content of the 
bracket is non-zero: 



Mi-* 2 )^ = ic, 



(18) 



where \C is the integration constant. Note also that, 
contrary to tpo, the sought for solution ipi is odd in x. 



Using the explicit form of ipo from Eq. ( |11| ) and inte- 
grating, one gets for x > 



tpx = CR 



dx = CR 



(1 



i dx. 



(19) 



r 



The integrals in expansion (|19|) are akin. Denoting 

F n = f x 2n e~ axl dx, 



one can easily write for them the recurrence relation and 
"initial" condition as 



F = ——F , 
da 



F n = 



vi(V^x), (20) 



respectively. Using the asymptotics of the error integral, 
with the exponential accuracy in a one finds 

F n = [(2n - 1)!!/2V"] F , F ~ ^/2^. (21) 



Comparing this with expression (|12|) for the function G, 
we get the representation 



ip 1 {x > 0) ~ CRF G. 



(22) 



Applying to Eq. (g2j) the normalizing condition ( 14 ) , one 
evaluates the constant as C — 1/RFqG. Therefore, from 



Eqs. (|20|)-(|22j) the principal rclaxational eigcnmode de- 
termined with the exp(— a) accuracy emerges as an odd 
step function 



¥>i(z) 



-1 for x < 0, 
1 for x > 0. 



(23) 



In Fig. [I] the limiting contour (^3|) is shown against 
the exact curves fi(x) obtained by solving numerically 
Eq. (||) for several values of a. We remark that in the 
statistical calculations carried out below the typical inte- 
grals are of two kinds. In the first, the integrand consists 
of the product of tpirfo and some non-exponential func- 
tion. As oc expert 2 , the details of behavior of (pi in 
the vicinity of x — are irrelevant because the approxi- 
mate integral will differ but exponentially from the exact 
result. The integrals of the second type contain dipi/dx 
in the integrand. For them a step-wise approximation 
Eq. (53) with its derivative equal identical zero every- 
where except for x = is an inadmissible choice. So, to 
keep the exponential accuracy in this case, one has to get 




FIG. 1: Eigenmode ipi(x) determined with the aid of the 
numerical solution of Eq. (|J) for the dimensionless barrier 
height a: 5 (dashed line), 10, 20, 25 (solid lines); the arrow 
shows the direction of a growth. Thick dashes show the step- 
wise function that is the limiting contour for ifi at a — » oo. 



back to Eq. (H|). 

The eigenvalue Ai corresponding to the approximate 
eigenfunction ipi from Eq. (|2|) is evaluated via for- 
mula (15) that can be rewritten as 



Ai = 



ipo 



2 1 



dx. 

(24) 



Substituting the derivative from Eq. (|18|), one finds 

Ai = C = (2/V5r) CT 1/2 /i?G, 
and using expression ( |T2] ) for R finally arrives at 
Ai = (4/xAF) a 3 / 2 e- CT /G 2 = A B /G 2 . 



(25) 



With G expanded in powers of cr _1 , see Eq. ([12|), this 
formula reproduces the asymptotic expression derived 
by Brown in Ref. |i~5l At G = 1 it reduces to his ini- 
tial result 



^ corresponding to the above-given Eq. (10). 
Function Ai (a) from Eq. (g5|) is shown in Fig. |^ in com- 
parison with the exact result obtained by a numerical 
solution. Indeed, at a > 3 the results virtually coincide. 

According to expansion (||), each decrement Ag defines 
the reference relaxation time 



Thence from Eq. (f2 



Tt = 2r D /X e 
we get 



t 1 = 2t d /\ 1 =t b G 2 



r B = 2r D /A E 



(26) 



(27) 



where tb denotes the asymptotic relaxation time ob- 
tained by Brown in Ref. [| Substituting in Eq. @ the 
explicit asymptotic series ( p"2| ) for G, one gets 



2a 3 / 2 



1 

a 



7 

4^2 



9 
2^3 



(28) 



FIG. 2: Asymptotic expression (|25|) for the eigenvalue Ai with 
allowance for terms up to <r~ 9 (solid line) compared to the 
exact numeric value (dashed line). 



D. Asymptotic integral time 

The decrements Ag or, equivalently, relaxation times Tg, 
being the characteristics of the eigenfunctions of the dis- 
tribution function, are not observable if taken as separate 
quantities. However, in combination they are involved in 
a useful directly measurable quantity, the so-called inte- 
gral relaxation time. In terms of correlation functions 
this characteristics is denned as 



Tint 



[m(t)m(0) ) 
(m 2 (0))o 



dt 



>(*MQ))o 

<z 2 (0))o 



dt. 



(29) 

where the angular brackets stand for the statistical en- 
semble averaging over the equilibrium distribution (p^). 
As follows from Eq. (|29|), the integral relaxation time 
equals the area under the normalized decay of magneti- 
zation. 

The Green function of Eq. (||), i.e., the probability den- 
sity of a state (x,t), provided the initial state is (xo,0), 
writes 



W(x,t;x ,0) = ijji(x) tpt(xo)e~ 



(30) 



1=0 



Similarly to Eq. (^), we expand the eigenfunctions in Leg- 
endre polynomials as 



^ = 1E(2H1)^( 



and introduce special notations for the first two functions 



k=l 



(32) 



fe=0 



rl>i = h J2(2k + l)Q k P k (x). 



k=0 



The procedures to evaluate the coefficients Sk and Q k 
and the explicit asymptotic forms for Qi and S2 are given 
in Appendix note representation ( |lT| ) for the equilib- 
rium function ipo. 
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Due to Eq. (|l4|), the coefficients in formulas ( |3l"| ) are terms one gets for the correlator in Eq. 



related to each other by b 



(PkPk' 



) a k , . 



In those 



((x(t)x(0))) o 



/ / xxp -0Q W(x, t; Xp, 0) dxdxp = \b^p 
J J l=1 



-\ t t/2T 



(33) 



r 



where the averaging over the current coordinate x is per- 
formed with the function W from Eq. ( |30| ) whereas that 
over the initial conditions — with the equilibrium function 
tJjq. Substituting expression ( |3^ ) in Eq. ( |29| ) one gets the 
integral time in the form 



Tint 



OO o OO OO n 



(34) 



Unlike n, which in principle cannot be evaluated) 18 
analytically at arbitrary a, for Tj n t an exact solution is 
possible for arbitrary values of the anisotropy parame- 
ter. Recently two ways were proposed to obtain quadra- 
ture formulas for T; nt . One method Jl9|| implies a di- 
rect integration of the Fokker-Planck equation. Another 
method p0[ involves solving three-term recurrence rela- 
tions for the statistical moments of W . The emerging 
solution for Tj nt can be expressed in a finite form in terms 
of hypergeometric (Rummer's) functions. Equivalence of 
both approaches was proven in Ref. EH]. 

In the present study, as mentioned, we are dealing in 
the high-barrier approximation. In this limiting case Ai 
is exponentially small, so that the term with I = 1 in 
the numerator in Eq. (p4]is far greater than the others. 
With allowance for Eq. (p2[) it can be written as 



/{'■ 



TiQl/(: 



(35) 



The equilibrium moment calculated by definition writes 
(x 2 ) Q = {l/2a){e a - 1) = l/G-l/2ff, (36) 
and for a » 1, using formula (AS) of Appendix |X] we get 



1/G. 



Substitution of Eqs. (p9) and (I 
for relationships (|l2|) , (p5|) and 
representation in the form 



(37) 

m i^[35|) with allowance 
27J) gives the asymptotic 




2oG 



Tint 



TB 



T~D 



{2a - G) 
^e° 



(38) 



1 



1 

a 



3 
2^2 



13 
4^3 



As it is seen from formulas ([28|) and ( |38| ) written with 
the accuracy up to cr -3 , the asymptotic expressions for 



the interwell and integral times deviate beginning with 
the term oc er~ 2 . This contradicts the only known to us 
asymptotic expansion of Tint given in Eq. (60) of Refs. 
and repeated in Eq. (7.4.3.22) of the book The 
latter expression written with the accuracy oc er~ 2 , in- 
stead of turning into Eq. ([38j) coincides with the Brown's 
expression (H|) for t\. Meanwhile, as it follows from for- 
mula ( J35| ) , such a coincidence is impossible and therefore 
Eqs. (60) of Refs. and (7.4.3.22) of Ref. || are mislead- 
ing. The necessity to rectify this issue made us to begin 
the demonstration of our approach with the case of the 
integral relaxation time. Further on we consistently ap- 
ply our procedure for description of the nonlinear (third- 
and fifth-harmonic) dynamic susceptibilities of a solid su- 
perparamagnetic dispersion. 



III. PERTURB ATIVE EXPANSIONS FOR THE 
DISTRIBUTION FUNCTION 

A. Static probing field 

To find the nonlinear susceptibilities, one has to take 
into account the changes that the probing field induces in 
the basic state of the system. In the limit a ^S> 1, which 
we deal in, the relaxation time T\ of the interwell mode 
ipi is far greater than all the other relaxation times Tfe. 
This means that with respect to the intrawell modes the 
distribution function is in equilibrium. So it suffices to 
determine the effect of the probing field H = Hh just on 
rpo and Assuming the energy function in the form 



U + U H = -Kv(e ■ nf - IvH{e ■ h), 



(39) 



[compare with Eq. @)], and separating variables in 
Eq. (3), one arrives at the eigenfunction problem 



Lf = ZVf , 



(40) 



where £ = IvH/kT and notation fp refers to the distribu- 
tion function modes that stem from ipo or ipi at H =/= 0, 
i.e., (3 = or 1. In Eq. ( [h)| ) operator L is defined by 
Eq. (§) whilst V = -£J(e x h) is the operator caused 
by the energy term Uh in (Eq. p%. As in above, for the 
non-self-conjugated spectral problem (|i(]) we introduce 
the family of conjugated functions gp and set fp = gpipo- 
Following our approach, in the low-temperature limit 
[a 1) we set to zero the eigenvalues corresponding to 
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both f and fx; compare with Eqs. ( O ) and Jl§| ) for 
and "01 . Assuming the temperature-scaled magnetic field 
£ to be small, we treat Uh as a perturbation Hamiltonian 
and expand the principal eigenfunctions as 

/o = E£"if\ /i = E^ l A (n) . (4i) 

71=0 71=0 

Thence for the field-free (H = 0) case one has = ipo 
and f[ 0) = ^i. The same kind of expansion is assumed 

for gp with g^ — 1 and g^ — tpx- Note also that in 
order to retain the normalizing condition we require that 
fp have zero averages. 

Substituting expansion (f4l| ) in Eq. (^o|) and collecting 

I 



the terms of the same order in £, we arrive at the recur- 
rence relation 

tff=Vff- l \ (42) 

that for the particular cases (3 = and 1 with the aid of 
the identity e x h = J(e • h) takes the formes 

J^ Q Jg { Q n) = Ji> g^ l) J{e ■ h), (43) 
Ji> Jgt l) = JtPog[ n ~ 1} J(e-h), 

respectively. Set (|||) solves easily for g since g^ = 
ipo = 1. Starting with n = 0, one gets sequentially 



o (1) 
yo 


= (e-h), 










(2) 
% 


= | [(e*) 2 " 


((e-/i) 2 )o] 








(3) 

9o 


= ^(e-^) 3 - 


i(e-h)<(e 








(4) 

So 


= M^) 4 


-((e-^) 4 > 




h) 2 ((e 


./i) 2 ) -((e.fc) 2 ) 2 ], 


(5) 

5o 


= IM( e ^) 5 - 




{{e-h?)o 




h) [((e./i) 4 ) -6((e 



(44) 



All the obtained functions are constructed in such a way that the corresponding /i satisfy the above-mentioned zero 
average requirement. We remark also that there is no problem to continue the calculational procedure to any order. 

Evaluation of g\ is done in two steps. At the first one, we set g^ equal to the antisymmetric step-wise function ( p^ ) 

and its derivative equal zero. After that from the second of Eqs. (^3|) we can express g[ k ^ in closed form. Taken up 
to the fourth order these "zero-derivative" solutions write 



9? 

(3) 
01 

(4) 

9x 



= tp x {e-h)- (<pi{e ■ h)) , 

= \tp\{e- Kf - (e- h) { <pi(e- h)) , 

= | [ Vl (e -hf-{ tpx (e . fc)3 )„] _ i( ^ ( C . h) > [(e • hf - ( (e ■ fc) 2 ) ] , 

= -^(e-fc) 4 - i(^(e-fc)> [(e./i) 3 -3(e-/i)(( e -/i) 2 )o] 

-i(e.^)<^ (e./i) 3 ; 



(45) 



Note the alternating parity in e with the term order growth in both Eqs. (|44|) and ([45| 



It is instructive to compare the approximate expres- 
sions ( f45| ) with the numerical results obtained without 
simplification of gf\ To be specific, we consider the 
case when probing field is applied along the particle easy 
axis n. Then Eqs. (^) become one-dimensional and the 
second of them writes 



(n) 



dx 



(n- 
9l 



(46) 



Its "zero-derivative" solutions up to the second order fol- 
low from the first two lines of Eqs. (Ma): 



9i 



(i) 



ipix- {<pix) , 



(2) 

9i 



\ifl x 2 - x( ifxx) ■ 



(47) 



In Figs. H and [| these functions are compared to the nu- 
merical solutions of Eq. (|4^). For our calculation, the 
most important is the behavior of those functions near 
x = ±1 since these regions yield the main contribution 
when integrated with the weight function ipg. As one 
can see from the figures, the "zero-derivative" solution 
g^ agrees well with the exact one, whilst g^ deviates 
significantly. This discrepancy is due to the change of 
the barrier height that occurs in the second order with 
respect to the probing field amplitude, and manifests it- 
self in all the even orders of the perturbation expansion. 
Correction of solution (|^J) makes the second step of our 
procedure. For that we integrate Eq. ( [46] ) two times by 
parts and substitute there the "zero-derivative" form of 
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0.4 



SS.2 



-0.5 



-0.2 



-0.4 



(2) 



0.5 1 



(2) 

FIG. 4: Function g\ found numerically (solid) and evaluated 
in the "zero-derivative" approximation (dashed). Asterisks 
show a corrected calculation with allowance for the coefficient 
L> 2 , see Eq. @. 



gP from Eq. @: 

9f ) = \x^ 1 -x{x^) + \ J x^dx, (48) 

(2) 

Thus one finds that the corrected g\ differs from this 



term 



9i 



(2) _ 



with the amplitude 
D 2 



1 

2./,, 



2%, 

x — — ax. 

ax 



(49) 



(50) 



We remark that the results of evaluation of the integrals 
I 2 k = f x 2k (dtpi/dx) dx can be arranged in the table 



(51) 



k 





1 


2 


I 2k 


1 


l-G- 1 


1-(1 + 1/2ct) G~ l 



so that Eq. (pfj|) gives 



Do 



1 



G- 1 



1 1 



37 



'2 A' 



. 2~7 2) 2G Aa ' 4cr 2 ' 8cr 3 ' 16cr 4 
Function g\ corrected in such a way is shown in Fig. H 

by asterisks. It is seen that the corrected dependence 
with a fairly good accuracy follows the the numerically 
obtained curve. 

In a similar way one can prove that the corrected func- 
tion g^ has the form 



(4) 



9i ' = ^ <Pl ^ - ~ [( <fi x ) x 3 - 3x ( x 2 ) ] 

(2) 

where the corrected function g\ given by Eq. 
used and 



(53) 



is 



10aG 2 ~ 22aG + G + 12a 



1 



32a 2 16a 3 32a 4 64<r 5 



1 



48aG 2 
5 29 



(54) 



In the general case, when the direction of the probing 
field does not coincide with the particle anisotropy axis, 

(«) 



of Eq. (|47|) by adding a step- wise [alike that of Eq. (23)] the corrected functions g\ still can be written as 



9i = \<P\ ( e ' h ) 2 - ( e ' h) (ifi(e ■ h)) + D 2 fi, 

ffl 3) = i (e • h) 3 - ( n (e ■ h) 3 ) ] - |(^i(e • h) } [(e • fc) 2 - ( (e ■ fc) 2 } ] + A^, 
ffl 4) = (e ./»)*_ (e - h) ) [(e ■ fc) 3 - 3(e -h)((e- h) 2 ) ] 

-i(e • h) ( <pi (e • h) 3 ) + I? 2 . g ( 2) + . (55) 



But since Eqs. ([43]) cannot be reduced to a form like solutions taking into account the behavior of function tpi 
Eq. (ff6|), the correcting coefficients D 2 and D4 cannot around zero are built up as power series near x — 0; such 
be presented in a closed form. In this case the corrected a procedure for the coefficients D 2 and D4 is described 



9 



in Appendix 



B. Dynamic probing field 

To obtain the dynamic susceptibilities, one has to find 
the distribution function W in the oscillating probing 
field £exp(ia;i). For this situation the kinetic equation 
(0) takes the form 



2T D ^+i\w(t) =£Ve i " t W(t), (56) 



Now we expand the functions subjected to operator L 
with respect to the set {ipk} of its cigcnfunctions, see 
Eq. (§): 



E4 1} 



(1) 



E 



f (i) 
Jo 



i>i ; 

(60) 

here (<p\f) denotes functional scalar multiplication, i.e., 
the integral of the product ipf over all the orientations 
of e. Substitution of Eq. ([30j) in (p3S|), multiplication of it 
from the left by ip^ and integration, render the expansion 
coefficient as 



where the operators L and V have been introduced in 
above. Assuming that the exciting field amplitude is not 
too high, we expand the steady-state oscillatory solution 
of Eq. (|5q) in a power series with respect to £ : 



W(t) = fW (n) e 4 ' 



(57) 



n=0 



Note that, mathematically, representation (pT\ ) is not 
complete. Indeed, in a general case the exact amplitude 
of the nw-mode must contain, along with the contribution 
~ £ n , an infinite set of terms ~ £ n+2 , £ n+4 j e ^ c - How- 
ever, in a weak field limit £ < 1 the terms with higher 
powers are of minor importance so that the main con- 
tribution to the magnetization response signal filtered at 
the frequency nu> is pro portional to £™. 

Substituting Eq. ( |57|) in j5q ) we arrive at the recur- 
rence set 



(58) 



that we solve sequentially starting from n = 1. At the 
first step the function in the right-hand side corresponds 
to the equilibrium case (£ = 0). Therefore, = tp , 

where the latter function is defined by Eq. ([ll]) and 
is frequency-independent. Combining Eq. (W2) written 
down for (3 = and n — 1 and Eq. (^8|), wc eliminate the 
operator V and get 



(2iujT D + L) W {1) = Lfk 



(i) 



(59) 



-fe M = (^|/o 1} ) I 1 + iuT k 
reference relaxation times 



(61) 
defined by 



where the reference relaxation times are 
Eq. @. 

In the low-frequency limit only cjti is set to be non- 
zero whilst all the higher modes are taken at equilibrium 
(cjTfe = 0). Thence, when constructing via Eq. d6C|), 

by adding and subtracting a term with c^(0), one can 
present the first-order solution in the form 



(i) 



IbJTi 
1 + lUJTl 



(^i|/o (1) ) 



Ipl 



(62) 



where /q , as seen from Eq. (|59|), is the equilibrium so- 
lution for the same value of the field amplitude £. We 
remind that the functions without upper index belong to 
the fundamental set defined by Eqs. (||) whereas those 
with an upper index are evaluated in the framework of 
the perturbation scheme described in Sec. Ill A. 



In the next order in £ the function is substituted 
in the right-hand side of Eq. ( p8| ) and through a proce- 
dure alike to that leading to Eqs. (|59|) - (|6l|) , the function 
i s found. We carry on this cycle up to k = 5. The 
results write 



ILUTl 
1 + iU)T\ 



(63) 



= /< 3) + (^|/} 2) ) - (H/o (3) )] /<°> - /< 

1 ^^i/^/^-i^i/^^i/fO/f 5 



(2) 



1 



4" iWTi 

1 

1 + 3zcijri 



(^|/o (3) ) + l(^|/^)(^lA (2) )]^ 



f(0) 



(64) 
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w& = C + 



+ IWTi 

1 

1 + SiLUTl 



^|/o (1) )a (33 -1(^|/o (1) )(^Ia (2) )a (1) 

(^|/o (3) ) + |(^|/^)(^IA (23 )]^ 



U)\ f(3) 



(1) 



(65) 



/ (5) -(^|/o (5) )a (0) + (^I/o (1) ) [(Ha (4) )a (0) -a ( 

Oi/o (1) ) (ha (2 0-(ha (3) )] [a (2) -(ha (2) )a ( 



(4) 
(0) 



1 



1 + ILOTi 
1 

1 + SiLUTi 

1 

1 + 5lU!Tl 



<Pi\fo 



f (4) 3 ( m I f (2)\ ,(2) 
A ~2 1^1 1 A J A 



A 



(0) 



¥ (^lA (2) ) 2 -t(^lA (4) ) 
(H/o (3) ) + 1 (^|/o (1) ) (H A (2) )] [a (2) - 1 (H A (2) ) A (0 

(^|/f) + i(^|/o (1 0(^lA (4 + l(^|/o (1 0(^lA (2) ) 

+§ (H/f ) (Ha (2) )] a ( 



(0) 



(66) 



We remark an important feature of Eqs. (J63|) — (J66|) : they 
do not contain dispersion factors of even orders. This 
ensures that the frequency dependence of the full distri- 
bution function W incorporates only dispersion factors 
with odd multiples of the basic frequency. Qualitatively, 
this is the result of absence of the interwell mode for 
the statistical moments of even orders. Technically, it 

is due to vanishing of the products yP\\f^\ entering 

Eqs. (|62|)-(|66|) if the sum k + I is even. This rule follows 
immediately from combination of the oddity of (fx, see 
Sec. U with the parity properties of the functions fjf^ 
introduced in Sec. Ill A. 



in the direction of the probing field. With representa- 



For actual calculations one needs the values of the 
scalar products entering Eqs. (|6^)-(|66|). In Appendix ^ 
we obtain their representations in terms of the moments 
Qk and Sk of the functions ipo and tpi j respectively. The 
procedures of asymptotic expansion of Qk and Sk are 
given in Appendix [A]. 



tion (57) for the distribution function, this magnetization 



component takes the form 



M = clv((e ■ h)) 



Tn+1 n+1 

71=1 



(kTy 



J (e h) W {n] de, 



(68) 



and the susceptibilities can be found by a direct compar- 
ison with Eq. (|67]). In other words, the set of x^ is ex ~ 
pressed through the perturbation functions found 
in the preceding section. Therefore, evaluation of x^ 
becomes, although tedious, but simple procedure. Re- 
markably, the final expressions come out rather compact. 



A. Linear susceptibility 



IV. DYNAMIC SUSCEPTIBILITIES 



The resulting expression can be presented in the form 



The set of magnetic susceptibilities of an assembly of 
non-interacting particles with the number density c is 
defined by the relation 



M 



X 



^H + x {3) H 3 + X {5) H 5 



(67) 



that describes the magnetization of the system in the di- 
rection of the probing field H = Hh. Therefore, of all the 
components of the corresponding susceptibility tensors, 
we retain the combinations that determine the response 



B, 



(i) 



D 



(i) 



1 + iuiT\ 



c/V 
3kT ' 



(69) 

which follows from substituting Eq. (|6^) in (|6q). Each 
of the two frequency-independent coefficients , being 
the result of statistical averaging over the orientational 
variable e, see Appendix ^|, expands into a series of Leg- 
endre polynomials with respect to /?, the angle between 
the direction h of the probing field and the particle easy 
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axis n. This can be written as 



=b^+b$P 2 (cosp), 
BW=&W + &Wp 2 (co S /3), 



(70) 



Definitions of functions Q\ and £2 and their explicit 
asymptotic representations are given in Appendix |a|. The 
asymptotic series for the coefficients derived on the 
base of expansion (021) and Eq. (B7h are 



6 (1) 

a 00 u 02 

6 (1) 

a 10 "12 



-0? 



2S 2 - 
2Q? 



1 _^ -.^ 34577 581133 

6 02 - „ + ,,„3 + Q-4 + 1R„5 + QO„6 + g 4cr 7 + 128(r 8 + " ' 



°10 



1 1 


13 


165 


2273 


- + TT + 
cr 4cr 3 




16<7 5 + 


32a 6 


1 3 


2 


31 


153 


a 4er 2 




~ 4^ ~ 


4ct 5 



3629 1564 785931 
16cr 6 ~ 64cr 8 



The other components, namely, byy and 6 12 ; , may be 
constructed straightforwardly using their relations with 
the given ones, see Eqs. (|70j). For a random system, 
that is for an assembly of non-interacting particles with 
a chaotic distribution of the anisotropy axes, the average 
of any Legendre polynomial is zero, so that BW^ = b^ , 
and the linear dynamic susceptibility reduces to 



1 + iujT\h 



'I Ml 



1 + jWTl 



(72) 



that is the asymptotic representation of the full expres- 
sion given by formula (39) of Ref. ||. 



triple frequency that at weak H scales as H 3 . Performing 
calculations along the same scheme as for x , one arrives 
at the sum-of-relaxators representation 



Y (3) _ 1 (3) [ B (3) 



(3) 



(3) 
X = 



d 4 v 4 



where the coefficients expand as 



B 



(3) 



1 + iu)T\ 1 + 'iiujT\ 



,(73) 



B. Cubic susceptibility 

As follows from definitions (^7|) and (|6^), the third- 
order susceptibility is defined through the response at the 



B k ] = b w + b k2 p 2 (cos p) + bfl Pi (cos 0), k = 0, 1 , 3 



up to the fourth Legendre polynomial in cos /?. 



(74) 



(3) 

The explicit expansions for the amplitudes b a L are 



1 47 49 815 7837 355391 

30a- 3 + 240a 4 + 40a 5 + 96a 6 + 120a 7 + 640a 8 + " ' ( ' 

1 2 4 1385 11231 19083 

42a 3 + 21a 4 + 7ct=" + 336a 6 + 336cr 7 + 64a 8 + ' ' ' 

2 8 41 50 1756 63749 



35a 3 35a 4 35a 5 7a 6 35a 7 160a 8 
1 1 23 61 1357 235447 11962691 694849241 15133953221 



+ . . . 



15 6a 240a 2 192a 3 960a 4 30720a 5 245760a 6 1966080a 7 5242880a 8 
13 65 25 863 3931 698911 35309123 2061480665 45071465669 



84 168a 168a 2 1344a 3 1344a 4 43008a 5 344064a 6 2752512a 7 7340032a 8 
112 1 73 17033 1007549 64390439 4493994417 



35 14a 35a 2 112a 3 560a 4 17920a 5 143360a 6 1146880a 7 9175040a 8 



+ . 



2 3 1 337 499 85309 2563751 245269747 47628510799 

15 + 10a + 16a 2 + 960a- 3 + 320a 4 + 10240a 5 + 49152a 6 + 655360a 7 + 15728640a 8 + " ' 
29 43 11 1279 1881 320765 48133699 920146163 178560431695 

84 + 56a + 56a 2 + 1344a 3 + 448a- 4 + 14336a 5 + 34406a 6 + 91750a 7 + 22020096a 8 + ' 

11 47 2 559 2419 409499 4080395 1166954357 75334335763 

105 + 210a + 21a 2 + 1680a 3 + 1680a 4 + 53760a 5 + 86016a 6 + 3440640a 7 + 27525120a 8 
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FIG. 5: Real (a) and imaginary (b) components of the cubic susceptibility of a superparamagnetic assembly with coherently 
aligned easy axes; the direction of the probing field is tilted with respect to the alignment axis at cos j3 — 0.5; the dimensionless 
frequency is ujtq — 10 -6 . Solid lines show the proposed asymptotic formulas taken with the accuracy cr~\ circles present the 
result of numerically-exact evaluation, dashed lin es co rrespond to the "zero-derivative" approximation (f45J). The discrepancy 
of the curves is commented in the text after Eq. (BIS). 



For a random system, the averages of Legendre polynomials drop out and BY 



,(3) 
J k0 



With respect to formalism 



constructed in Ref. H, the above expressions yield the asymptotic representations for formulas (42) and (43) there. 



C. Fifth-order susceptibility 



The susceptibility of the fifth order writes in an expectable way as a sum of three relaxators: 



y (5) _ ± (5) B (5) -1 , _ 



B 



(5) 



B 



(5) 



B 



(5) 



+ Huit 1 + hlUJT 



(5) CI 6 V 6 

Xo - 



(fcT) £ 



(76) 



with the coefficients 



B k ] = b m + b k2 P 2(cos(3) + 6j^P 4 (cosj8) + b^P 6 (cos0) k = 0, 1,3, 5. 



(77) 
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The explicit asymptotic series are 



367 123 41233 



2240a 6 70a 7 2240a 8 
19 1 65 79 85913 72636131 4543038053 14938598691 



96 


420a 


120a 2 


47 


11 


29 


~560 


+ 35a 


280a 2 


311 


137 


13 


3360 


420a 


+ 105<r 2 



+ . . . 

1920a 3 ' 4480a 4 ' 143360a 5 1 3440640a 6 ' 1835008a 7 ' 73400320cr 8 + ' ' 
5911 2141 1874309 299470403 17964831133 400677748549 



1792<t 3 336cr 4 57344cr 5 6881280a 6 55050240cr 7 20971520a 

437 5473 1046209 169435283 684614895 230861266333 



26880a 3 1920<T 4 286720cr 5 6881280a 6 55050240a 7 146800640cr 8 



507 5377 
+ ■ 



112cr 5 28a 6 448ct 7 448a 8 
13 19 23 737 2959 99733 50499149 350973527 72765921299 

504 168ct + 672a 2 8064a 3 5376a 4 28672a 5 2064384a 6 1835008a 7 44040192a 8 + " ' 
5 149 193 5245 18677 1785635 289305193 5846947361 394448762615 

21 168ct 672a 2 8064a 3 5376cr 4 86016a 5 2064384a 6 5505024a 7 44040192a 8 
139 27 109 1343 9203 431321 9839105 196654913 30690812563 

504 " 28a + 336a 2 2016<r 3 ^ 2688a 4 _ 21504a 5 _ 73728a 6 ~ 196608a 7 3670016<t 8 + " ' 

3 1563 7767 613353 

~ 140a 5 ~ 6160a 6 ~ 3080<J 7 ~~ 24640a 8 + " ' 
15 183 713 433 409 319665 8222083 5744848239 403943151013 



2464 6160a 24640a- 2 19712cr 3 4928c 4 630784a 5 2293760<t 6 201850880a 7 1614807040a 8 

29 293 47 7081 74647 7137293 385804437 4682760003 1580817298041 

~ 280 + 770a 385a 2 + 24640a 3 + 49280a 4 + 788480a 5 + 6307840a 8 + 10092544a 7 + 403701760a 8 + ' " 

1713 2929 2551 34863 92061 34432191 23756287 15647080587 7317549380671 



12320 6160a 24640a 2 98560a 3 49280a 4 3153920a 5 327680a 6 28835840a 7 1614807040a 8 
1 3 1467 



616a 6 77a 7 2464a 8 

11 7 337 53 51433 2188103 471762913 4428495037 

~ 1584 + 1848a + 1056a 2 + 88704a 3 + 1848a 4 + 315392a 5 + 2064384a 6 + 60555264a 7 + 69206016a 8 + ' " 
1 10 17 103 2489 236615 38344237 776232845 52467158027 

84 + 231a 924a 2 + 3168a 3 + 14784a 4 + 236544a 5 + 5677056a 6 + 15138816a 7 + 121110528a 8 + ' " 
1207 661 17 5525 1169 9116467 131486063 1918435847 639291980689 

55440 9240a 12320a 2 88704a 3 3520a 4 4730880a 5 10321920a 6 20185088a 7 807403520a 8 



(78) 



and for a random system, as for the lower orders, B k = 



V. DISCUSSION 

The above derived formulas despite their hefty look 
are very practical. Indeed, they present the nonlinear 
initial susceptibilities of a superparamagnetic particulate 
medium as analytical expressions of arbitrary accuracy. 
With respect to the frequency dependence they give the 
exact full structure of the susceptibility and prove that 
it is very simple thus putting former intuitive considera- 
tions on a solid ground. This makes our formulas a handy 
tool for asymptotic analysis. Yet more convenient they 
are for numerical work because with their use the difficult 
and time-consuming procedure of solving the differential 
equations is replaced by a plain summation of certain 
power series. For example, if to employ Eqs. (|72|)-(f78|), 
a computer code that fits simultaneously experimental 
data on linear and a set of nonlinear susceptibilities tak- 
ing into account the particle polydispersity of any kind 
(easy axes directions, activation volume, anisotropy con- 
stants) becomes a very fast procedure. 

Graphic examples justifying our claims are presented 



in Figs. g| and |6J, where the components of two nonlinear 
complex susceptibilities are plotted as the functions of 
the parameter a. For a given sample, a in a natural way 
serves as a dimensionless inverse temperature. In those 
figures, the solid lines correspond to the above-proposed 
asymptotic formulas where we keep the terms up to <r~ 3 . 
The circles show the results of numerically-exact solu- 
tions obtained by the method described in Ref. || Note 
that even at a ~ 5 the accuracy is still rather high. 

The model that may be called the predecessor of the 
afore-derived results was proposed in Ref. ^3|. There, 
the authors calculated the initial susceptibilities up to 
the seventh order having replaced a superparamagnetic 
assembly by a two-level macrospin system. The inter- 
relation between the present work and Ref. ^3] closely 
resembles the situation with the evaluation of the rate of 
a superparamagnetic process. First in 1949 Neel[Q and 
then, ten years later, Brownj3j had evaluated the super- 
paramagnetic time in the framework of a two-level model. 
In such a framework, one allows for the magnetic moment 
flips but totally neglects its possible diffusion over ener- 
getically less-favorable directions. In 1963 BrownH had 
succeeded to overcome this artificial assumption and took 
into account the possibility for the magnetic moment to 
wander over all 47r radians. 

In the present case, the obtained v/T dependencies of 
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given in an analytic form. The authors propose to eval- 
uate them by solving an infinite set of recurrence equa- 
tions. Hence, the procedure! 23 does not provide any gain 
with respect to former ones neither in analytical consid- 
erations nor in constructing fitting codes. 

In the presented framework the results by Klik and 
Yao (including the analytical formulas for them missing 
in Ref. ||^) can be obtained immediately if to take the 
function (p% in a step- wise form ( p3| ) and not to allow for 
the corrections caused by the finiteness of its derivative at 
x = 0. In our terms this means to stop at set (|45|), i.e., 
"zero-derivative" solution, and not to go further. The 
emerging error is however, uncontrollable and not at all 
small. As an illustration, in Fig. |B| we show the result 
obtained with this model (dashed lines) for the cubic sus- 

(3) 

ceptibility X3J m a textured system where the particle 
common axis n is tilted under the angle (3 = 7r/3 to the 
probing field. One can see that deviations are substan- 
tial. 

In Ref. | we have proposed, although without rigorous 
justification, a formula for the cubic susceptibility of a 
random assembly 
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(3) _ 



1 (3) (l + 25 2 2 )(l-^r 1 ) 



45(1 + zwri)(l + 3?wti) 



(79) 



FIG. 6: Real (a) and imaginary (b) components of the fifth- 
order susceptibility of a random superparamagnetic assembly; 
the dimensionless frequency is u)To = 10~ 6 . Solid lines show 
the proposed asymptotic formulas with the accuracy er~ 3 , cir- 
cles present the result of a numerical evaluation 



the nonlinear susceptibilities and those from Ref. [2^ are 
qualitatively the same. Their most typical feature is the 
double-peak shape. Quantitatively, however, the corre- 
sponding lines differ and do not reduce to one another 
in any case. Indeed, as long as the temperature is finite 
(whatever low) , the configurational space for the unit vec- 
tor e of the magnetic moment is the full (47r-radian) solid 
angle; its reduction to just two directions along a bidirec- 
tional axis could not be done otherwise than "by hand" . 
This is exactly what the two-level Ising-like model does: 
it forcibly imparts a quantum property (discrete spin pro- 
jections) to a macrospin assembly. From the calculational 
viewpoint, another essential demerit of the results |23j is 
that the coefficients in the susceptibility formulas are not 



that proved to be well adjusted for approximating the 
results of numerical calculations in all the temperature 
interval and also appeared to be good for fitting exper- 
imental data.Q Now we see that this very expression 
follows from Eqs. fl73|)-(f75|) if to expand the coefficients 

b^jn up to the zeroth order with respect to a^ 1 . This jus- 
tifies Eq. ((79|) as a formula yielding a correct frequency 
dispersion of the cubic susceptibility of a random assem- 
bly at low temperatures. The cause of its applicability 
at high temperatures is the exponential dependence of T\ 
on a. Indeed, in the frequencies range cotq -C 1, where 
we work, the condition a 25 1 means t± — > tq, and all 
the dispersion factors in Eq. d79) drop out. This trans- 



forms expression (79) in a correct static susceptibility 
that is also a true result. To avoid any confusion we re- 
mark that Eq. (|79| ) differs from formula for \^ given in 
Ref. by the coefficient (—1/45) due to the difference 

in definitions: in Ref. Q it was included in \^ . 

Applying the similar procedure to Eqs. (f7q)-(f78|) we 
get the expression for the fifth-order susceptibility 



(5) 



1 (5) (2 + 12gf + 4gf) 



21 



16 



X6 



ILUTl 



945 



(1 + iu!Tl)(l + 3iU!Tl)(l + 5iu>Ti) ' 



(80) 



that, following the example of the already tested Eq. (79), in the whole temperature interval. As we have already 
has high chances to be a good approximation for x£% ascertained in Ref. §, the best interpolation expression 
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for the relaxation time in the susceptibility formulas is 

e CT - 1 f cr fa „ ,1 _1 
2er [ 1 + cr V 7r 

proposed in Refs. [||[ |l4| . 



The package Maple V for PC used in the calculational 
work was obtained by Institute of Continuous Media Me- 
chanics in the framework of the EuroMath Network and 
Services for the New Independent States - Phase II (Em- 
Net/NIS/II) project funded by INTAS under grant No. 
IA-003. 



VI. CONCLUSIONS 

A consistent procedure yielding the integral relax- 
ation time and initial nonlinear susceptibilities for an as- 
sembly of non-interacting superparamagnetic particles is 
constructed in the low-to-moderate temperature range. 
Starting from the micromagnetic kinetic equation that 
describes intrinsic rotary diffusion of the particle mag- 
netic moment, we obtain the results in an analytical 
form. They are presented as asymptotic series with re- 
spect to the dimensionless parameter a that is the uni- 
axial anisotropy barrier height scaled with temperature. 
High-order expansion terms are easily accessible that al- 
lows to achieve any desirable extent of accuracy. This 
is proven by comparison of the proposed approximation 
with the numerically-exact results. The susceptibilities 
contain angular dependencies that allow one to consider 
the particle assemblies with any extent of orientational 
texture — from perfectly aligned to random. The new for- 
mulas stand closer to reality than those for a two-level 
system and are to facilitate considerably both analytical 
and numerical calculations in the theory of superparam- 
agnetic relaxation in single-domain particles. 
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APPENDIX A: EVALUATION OF THE 
EXPANSION COEFFICIENTS FOR 
EIGENFUNCTIONS Vo AND Vi 

Both functions ipo and ip\ are uniaxially symmetrical 
about the anisotropy axis n and can be expanded in the 
Legendre polynomial series, see Eq. (|32|): 

ip = lJ2{2k + l)S k P k (x), k = 0,2,4,..., (Al) 

fe=0 

V>i = I £ (2* + l)Q k Pk(x), k = 1, 3, 5, . . . , 
fe=i 

where in accordance with the parity properties of the 
eigenfunctions non-zero terms are 

S = l, S k = (P k (x)\tp ), ft = 2,4, . .., (A2) 

Qk = (Pk(x)\ipi) , k = 1,3,5,.... 

Taking into account that -^l = Vwii where tpa i n a fi- 
nite form is given by Eq. ( |ll| ) , one arrives at the general 
formula 

T k = (1/R) f P k {x)e ax2 dx, (A3) 
Jo 

where T is S k for even and is Q k for odd values of the 
index, and the function R(a) is defined by Eq. ([Tl|). In 
particular 

Q l = {l/R) xe™* dx=\{e a -l)/aR. (A4) 
Jo 

Using asymptotic expansion (|l^) for R, one gets 



1 1 5 37 353 4881 55205 854197 
Ch - 1/G - 1 - — - ^2-^3-^4 - - - - "128^8" + ■ ' ' ^ A5 ^ 

Knowing Q\, one can derive all the other moments Q k with the aid of the three-term recurrence relation obtained 
from Eq. (||) by setting there b k = Q k and A = 0. The same relation can be used to find the equilibrium order 
parameters S k . This is a head-to-tail procedure, where Sq = 1 and S2 is determined by the integral 

pi 

S 2 = (1/2R) / (3x 2 - 1) e ax dx. (A6) 



Taking the latter by parts one gets 



S2 = 1 [e a - R] /*R. 
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On comparison with Eq. (A4), we find 



5' 2 



(3 - 2a)/a, 



that upon substituting asymptotic series (A5), transforms into 



So = 1 - — - 



2<T 4cr 2 



15 

8^ 



111 

16(T 4 



1059 12243 165615 2562591 



32ct 5 

r 



64cr b 



128cr' 



256ct 8 



(A7) 



APPENDIX B: EVALUATION OF THE 
CORRECTING COEFFICIENTS D n IN A 
GENERAL CASE 



Let us present the solution of Eq. (42) in the form 

S n \ (Bl) 



where the functions are rendered by formulas J45| ) 
and are not corrected with respect to the derivative 
difx/dx. Substituting Eq. (Bl) in ( p2| ) and taking into 
account Eqs. (|4^), we get a recurrence sequence of equa- 
tions for the corrections vS n ': 



With allowance for the fact that function tp^ depends 
only on x, Eq. (B2) rewrites as 



(B2) 



dx 



^o(l - x 4 ) 



2 (e • h) n dtpi 



dx 



Finally, making use of the relation 

dpi Ai 
dx 2-00 (1 — x 2 

that follows from Eq. (Eq), we get 



2 dx 



(e ■ hy 



In particular, at n = 1 Eq. (B4) takes the form 



Ai d 



2 da; 



(e-h) 



(B4) 



(B5) 



Equ atio ns (B4) are solved sequentially beginning from 
Eq. (B5) by expanding in a power series with respect to 
x. The right-hand sides of Eqs. (B4) and (B5) are pro- 
portional to an exponentially small parameter Ai. Just 
due to that we did not take into account the corrections 
of the order tt(") when deriving Eqs. ([45]) ■ However, the 
quantities 



D n = 



n = 2,4.. 



have finite values. To show that, let us multiply Eq. (B4) 
by ipi and integrate. This yields 



(e • h) n 



(B6) 



In the left part we make use of the fact that tpi is the 
left eigcnfunction of the operator L, in the right part the 
integrals are taken by parts and yield 



AiA, 



2 / (1 



x 2 )u( n -V 



dtpi 



dx 



dx 



(e-h)dx (B7) 



1 dipx (e • h) r 
dx nl 



dx. 



Replacing the derivative d<p\/dx in the first term of the 
right-hand side with the aid of Eq. (B3), we arrive at the 
representation of the coefficient D n as 



D n = 



1 d 



ij)Q dx 



{e . h)dx - I 1 ^fl^L dx . (as) 
o dx n! 



Since ipo exp(crx 2 ), the first integral in Eq. (B8) can be 
presented as an asymptotic series if the power expansion 
of the function u^" 1 ) in the vicinity of x = is known. 
A closed form for the second integral can be found with 



(B3) the aid of the table given in Eq. (|51|), see section III A 



As an example, we calculate the coefficient D^. 
from the addition theorem 



(e • h) = cos 6 cos (3 + sin 8 sin f3 cos ip, 



As 



we seek the solution of Eq. (|B5| ) the sum 

u« = cos (3j2 C k^ k + sin/?e^(l - x 2 )^ C^x k . 

k k 

(B9) 

Here the upper index of the C coefficients corresponds to 
the azimuthal number m of the spherical harmonic e lmv . 
Operator L now includes the azimuthal coordinate and 
takes the form 



- L = (1 - x 2 )^ - [2ax(l - x 2 ) + 2x] — 



+ 



2a(3x 2 - 1) 



m 



l-x 2 



Substitution of expansion (B9) in Eq. (B5) leads to the 
set of equations 
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2a(k + m + 1)C£ 



(m) 
2 



[fc(fc + 1 + 2m + 2a) + m(m + 1) + 2a] C ( ™ ] + (Jfe + l)(jfc + 2)C££> 



iV. 



(m) 



(BIO) 



r 



where to = 0, 1 and the numbers in the right-hand side 
are 



N; 



(0) 



-1 for k = 0, 
for jfe ^ 0, 



iV, 



(i) 



1, for fe odd, 
0, for k even. 



In reality, one retains in expansion (B9) only a finite 
number of terms so that Eqs. ( B10| ) could be easily solved 
analytically by any computer algebra solver. In terms of 
expansion (JB9|) expression (JB8|) at n = 2 writes 



£ 2 = cos 2 /?£cf fe 



( 0) (2fc-l)!! 



2 k a k G 



(Bll) 



isin^^cW, 
fc=i 
1 2G- 3 
~6 6G - 



(1) (2fc-l)H 
2 k a k G 

P 2 (cos/3). 



Since the coefficients C found from Eq. ( B10 ) are func- 
tions of a, one has to perform in Eq. (Bll) asymptotic 
expansion. This gives finally 



1 



Do = — 



sin z P 



1 

4.rr ' 4<7 2 
1 1 



5 

8a 1 ' 
1 



37 



16a 4 



(B12) 



7 
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4 4 8a 16a 2 64a 3 64a 4 ' 
As it should be, at j3 = this formula reduces to Eq. 15 



that was obtained for a one-dimensional case. We re- 
mark, however, that in a tilted situation (/3 ^ 0) the 
coefficient D 2 acquires a contribution independent on a 
that assumes the leading role. This effect is clearly due 
to admixing of transverse modes to the set of eigenfunc- 
tions of the system, and it is just it that causes so a 
significant discrepancy between the "zero-derivative" ap- 
proximation and the correct asymptotic expansion for 
-^(3) curves i n pig H Evaluation of the coefficient D4 is 
done according to the same scheme and requires taking 
into account a number of the perturbation terms that 
makes it rather cumbersome. 



APPENDIX C: EVALUATION OF INTEGRALS 

Before proceeding to the integrals (scalar products) in 
Eqs. (|62])-(|66|) and (|6g|), let us consider the "primitive" 
ones 



((e-h) n \ip ), Y n = ((e-h) n \*p 1 ) 



The functions 4>o an d are originally defined in terms 
of the angle 9 = arccos(e ■ n). Thus, before performing 
integration one needs to transform both integrands to 
the same set of angles. Doing this with the aid of the 
addition theorem for Legendre polynomials, one finds 



and 



X 2 = § [2S* 2 P 2 (cos/3) + 1], X 4 = i [8S 4 P 4 (cos /3) + 20S* 2 F 2 (cos/3) + 7] 
X 6 = [16S 6 P e (cos P) + 72S 4 P A (cos (3) + 11QS 2 P 2 (cos P) + 33], 

Y 1 = Qi cos P , Y z = \[ 2Q 3 P 3 (cos P)+3Q 1 cos p] , 
Y 5 = i [8Q 5 F 5 (cos/3) + 28Q 3 P 3 (cos/3) + 27Q 1 cos/3], 



(CI) 



(C2) 



where cos/3 — (n ■ h) and the parameters Sk and Qk are the expansion coefficients introduced by formulas (f32|). 

Now using the expressions for functions /q 7 ^ and f[ n ^ derived in Sec. Ill A one sees that the relevant integrals of 
Eqs. (|63|)-(|68|) are expressed in terms of Xk and Yfc as 



((e • h)|/ (1) ) = X 2 , ((e ■ h)\f^) = \X 4 - iXl ((e • h)\f^) = ^X e - \X A X 2 + IXf; 



(C3) 



(fi\fi 1] ) = Yi , = $Y 3 - \X 2 Y X , (y x \ff) = j^Y 5 - ±Y 3 X 2 + \XlY x - ±X 4 Y t ; (C4) 
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(</>i|/i (2) ) = 5^2 - Y? + D 2 , ((e • fc)|/< 2 >) = ir 3 - Y,X 2 + D 2 Y l , 

(^l|/l (4) ) = 21^4 ~ |i3il + \-X 2 Y? + £>4 + D 2 (H/l (2) ) ' 

((e • h)|/ x (4) ) - ^y s - ir 3 X 2 - \Xtf x + \XlY x + D 4 Yi + Da ((e • fr)|jf ] ) . (C5) 
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